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A simple trajectory model has been developed and is presented. The particle trajectory path is 
estimated by computing the vertical position as a junction of the horizontal position using a constant 
horizontal velocity and a vertical acceleration approximated as a power law. The vertical particle 
position is then found by solving the differential equation of motion using a double integral of vertical 
acceleration divided by the square of the horizontal velocity, integrated over the horizontal position. 
The input parameters are: Xoand yo, the initial particle starting point; the derivative of the trajectory at 
Xoand yo, So - s(xo)= dx(y)/dy\ y=y o ; and b where bxo/yo Is the final trajectory angle before gravity pulls 
the particle down. The final parameter Vo is an approximation to a constant horizontal velocity. This 
model is time independent, providing vertical position x as a function of horizontal distance y: 
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The first term on the right in the above equation is due to simple ballistics and a spherically expanding 
gas so that the trajectory is a straight line intersecting (0,0), which is the point at the center of the gas 
impingement on the surface. The second term on the right is due to vertical acceleration, which may be 
positive or negative. The last term on the right is the gravity term, which for a particle with velocities 
less than escape velocity will eventually bring the particle back to the ground. The parameters b, So, 
and in some cases vq, are taken from an interpolation of similar parameters determined from a CFD 
simulation matrix, coupled with complete particle trajectory simulations. 
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Introduction 

In preparation for the Apollo program, Leonard Roberts of the NASA Langley Research Center 
developed a remarkable analytical theory that predicts the blowing of lunar soil and dust beneath a rocket 
exhaust plume. Roberts assumed that the erosion rate was determined by the excess shear stress in the gas 
(the amount of shear stress greater than what causes grains to roll). The acceleration of particles to their 
final velocity in the gas consumes a portion of the shear stress. The erosion rate continues to increase until 
the excess shear stress is exactly consumed, thus determining the erosion rate. Roberts calculated the 
largest and smallest particles that could be eroded based on forces at the particle scale, but the erosion rate 
equation assumed that only one particle size existed in the soil. He assumed that particle ejection angles 
were determined entirely by the shape of the terrain, which acts like a ballistic ramp, with the particle 
aerodynamics being negligible. The predicted erosion rate and the upper limit of particle size appeared to 
be within an order of magnitude of small-scale terrestrial experiments but could not be tested more 
quantitatively at the time. The lower limit of particle size and the predictions of ejection angle were not 
tested. 

We observed in the Apollo landing videos that the ejection angles of particles streaming out from 
individual craters were time-varying and correlated to the Lunar Module thrust, thus implying that 
particle aerodynamics dominate. We modified Roberts’ theory in two ways. First, we used ad hoc the 


ejection angles measured in the Apollo landing videos, in lieu of developing a more sophisticated method. 
Second, we integrated Roberts’ equations over the lunar-particle size distribution and obtained a compact 
expression that could be implemented in a numerical code. We also added a material damage model that 
predicts the number and size of divots which the impinging particles will cause in hardware surrounding 
the landing rocket. Then, we performed a long-range ballistic analysis for the ejected particulates. 

An ongoing activity since that time has been the development of particle trajectory simulations in order to 
improve estimates of trajectory angle and speed as a function of particle size and initial starting point on 
the ground relative to the engine nozzle. Four distinct models have been used to attack this problem: 

1. Computational Fluid Dynamics (CFD) code ( Fluent ™) - the simulation is set up to characterize 
the gas flow from a rocket engine impinging vertically on a rigid horizontal surface. The surface 
may have other specific features that replicate shallow circular craters or vertical berms to stop 
the high velocity spray of regolith from the surface. 

2. Particle Trajectory with Quadratic Shepard's Interpolation (PTQ) - a stepped integration method 
to solve a single particle’s differential equations of motion. The CFD computed gas properties, 
temperature, density, and vector velocity, are the input to this model, where the gas provides the 
force to propel a single particle of diameter D from the ground at an initial starting point, r 0 . 

3. Empirical Trajectory Path (ETP) model - an empirical time independent function is used to 
estimate a 2D trajectory path of a single particle, including path shape and angle. The horizontal 
velocity v y , initial starting point (jc 0 , yoX and two shape parameters b and c constitute the model’s 
input parameters.. The shape parameters are determined by fitting the empirical function to 
associated PTQ trajectories from step 2. A matrix of ETP shape parameters and horizontal 
velocities are generated from a matrix of particle variables, such as of particle size, starting 
locations, engine height, and engine thrust, and by fitting the ETP model to the PTQ generated 
trajectories. Using a simple A-dimensional interpolation algorithm, ETP model parameters can 
then be estimated for any particle size, starting point, as well as engine height above the ground. 

4. Physical Trajectory Path (PTP) model - The same overall procedure to obtain parameter 
interpolations used in step 3 is used in this case, but the empirical function is replaced with a 
special solution to the differential equations of motion. This solution is based on an estimate of 
vertical particle acceleration and particle horizontal velocity. This model uses two shape factors, 
So and b , as well as horizontal velocity v 0 and initial starting point (jc 0 , yo)- The only difference in 
the parameter definitions between ETP and PTP is the parameter c vesus s 0 . 


Empirical Model 

In previous work (Lane, 2010), an empirical trajectory path model was investigated to describe the 
trajectory path of particles on the Lunar surface under a rocket plume: 


X(y ) = bx 0 

with a gravity free ( g = 0) trajectory angle: 
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The model of Equation (1) mimics the general trajectory path behavior of particles on a smooth horizontal 
surface propelled by a vertically impinging rarefied gas jet. Figure 1 shows a plot of Equation (1) where x 
is the vertical direction (direction of gravity), y is the horizontal direction, and (jc 0 , yo) is the initial starting 
point of the particle. When b - 1, with gravity g = 0, the resulting path is that of a straight line which 
intersects the points (0, 0) and (x 0 , yo)- The trajectory angle in this case is tan _1 (xo/yo). In most cases the 
parameter b ^ 0 so that the final gravity free trajectory angle is tan ~\b x 0 /yo). The parameter c in Equation 
(1) adjusts the curvature ofx(y). 

The “gravity free” region in this context is determined by the last term on the right hand side of Equation 
(1). When this term is much smaller than the first term on the right hand side, the trajectory angle will be 
a constant value equal to tan _1 (£> x 0 /y 0 ). The trajectory path model parameter v y of Equation (1) is the 
horizontal velocity of the particle. When this value is very large, from 100 to 2400 (escape velocity) m/s, 
the gravity term (g = 1.622 m s" 2 on the lunar surface) has minimal effect within a radius of 20 to 30 m for 
smaller particles less than 100 |Ltm in diameter. However, regardless of the values of g or v y the straight 
trajectory path defined by tan _1 (Z? %o/y 0 ) will eventually lose out to the gravity term and the path will curve 
back to intersect the surface (unless the particle velocity exceeds escape velocity). 



Figure 1. Solid thick line: example trajectory path from Equation (1) with x 0 = 0.02 [m], y 0 = 1.0 [m], b = 
10, c = 1, v y = 10 [m/s]; solid line: gravity term (last term on right hand side of Equation (1); dashed line: 
tan _1 (xo/y 0 ) slope line; dotted line: tan _1 (Z? x^lyj) line; thick dashed line: gravity free version of Equation 

(1), with g = 0. 


CFD Trajectories Interpolation Matrix 

Table 1 shows a matrix of N particle values (size and starting location), P t = (D h x 0l , yoi) T , with the 
corresponding of fit of Equation (1) to the final CFD based trajectories. The fitted model parameters are 
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The values in Table 1 are found by first computing a matrix of full particle trajectories generated by post- 
processing CFD output data using the equations of motion, Equation (1.0) and considering the details of 
particle drag and lift (Lane, 2010). Figure 2 shows trajectories generated by fitting the model of Equation 
(1) CFD based trajectories. The result of the fit for N trajectories is shown in Table 1. 


Table 1. Empirical model parameters found by fitting Equation (1) to a 
matrix of CFD trajectory solutions for a range of particle sizes D and 
initial starting points x 0 , yo- 


D 

x 0 

Vo 

b 

c 

b 

1 

0.01 

0.88779 

5.061 

0.6445 

1997 

1 

0.01 

1.08699 

4.715 

0.7719 

1770 

1 

0.01 

1.33089 

-99 

-99 

-99 

1 

0.01 

1.62951 

-99 

-99 

-99 

1 

0.01 

1.99514 

8.221 

0.58 

640.4 

1 

0.01 

2.4428 

-99 

-99 

-99 

1 

0.01 

2.99091 

-99 

-99 

-99 

1 

0.01 

3.662 

24.78 

0.7374 

414.4 

1 

0.01 

4.48367 

312.7 

1.163 

480.3 

1 

0.01 

5.4897 

-99 

-99 

-99 

1 

0.01 

6.72146 

49.69 

1.269 

346.1 

1 

0.1 

0.88779 

1.953 

0.5507 

2381 

1 

0.1 

1.08699 

-99 

-99 

-99 

1 

0.1 

1.33089 

2.254 

0.7018 

2228 

1 

0.1 

1.62951 

-99 

-99 

-99 

1 

0.1 

1.99514 

2.402 

0.7626 

1872 

1 

0.1 

2.4428 

-99 

-99 

-99 

1 

0.1 

2.99091 

2.42 

0.8279 

1098 

1 

0.1 

3.662 

-99 

-99 

-99 

1 

0.1 

4.48367 

3.123 

0.7227 

494.7 

1 

0.1 

5.4897 

-99 

-99 

-99 

1 

0.1 

6.72146 

5.356 

0.8368 

311.8 

1000 

0.01 

0.88779 

3.177 

0.7301 

84.87 

1000 

0.01 

1.08699 

-99 

-99 

-99 

1000 

0.01 

1.33089 

-99 

-99 

-99 

1000 

0.01 

1.62951 

4.057 

0.8624 

41.21 

1000 

0.01 

1.99514 

4.296 

0.7009 

26.22 

1000 

0.01 

2.4428 

10.78 

0.5 

17.12 

1000 

0.01 

2.99091 

-99 

-99 

-99 

1000 

0.01 

3.662 

1.739 

2 

4.296 

1000 

0.01 

4.48367 

-99 

-99 

-99 

1000 

0.01 

5.4897 

-99 

-99 

-99 

1000 

0.01 

6.72146 

1 

0.5 

0 

1000 

0.1 

0.88779 

1.591 

0.4195 

133.6 

1000 

0.1 

1.08699 

-99 

-99 

-99 

1000 

0.1 

1.33089 

2.136 

0.7379 

116.3 

1000 

0.1 

1.62951 

-99 

-99 

-99 

1000 

0.1 

1.99514 

2.289 

0.8039 

88.73 

1000 

0.1 

2.4428 

-99 

-99 

-99 

1000 

0.1 

2.99091 

2.239 

0.8885 

43.09 


Vertical Distance from Surface [m] Vertical Distance from Surface [m] 




Figure 2. Solid lines: Full particle trajectories generated from post-processing of CFD simulations; dashed 
lines: trajectories generated by fitting the model of Equation (1) to solid lines. Top: D = 1 pm and x 0 = 

0.01 m; bottom: D = 1000 pm and xq = 0. 1 m. 


Note that an entry of “-99” in Table 1 corresponds to no fit data available. Each entry in Table 1 is the 
result of a full CFD based trajectory, but not all of those entries were fitted to Equation (1). The number 
of useable entries N from Table 1 is therefore 21, in this example. 


Once a matrix of model parameter values has been determined, an estimated set of parameters P can be 
computed for any arbitrary set of input values X using an application of Shepard’s interpolation algorithm 
(Shepard, 1968) to the matrix defined by Table 1: 


P = 
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where, 
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Theoretical Model 


The general equation of motion of a particle of mass m acted upon by a force F is: 

mr = F . (6) 

Since we are interested in the trajectory path P(y), the following substitution eliminates the time variable 
from the differential equation of Equation (6): 
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Note that in Equation (7), the independent variable is y while the dependent variable is x = P(y). This 
seemingly backwards convention comes about because the axis of the gas jet is defined to be along the 
vertical x-axis while the ground plane lies along the y-axis. Considering only the 2D case (the following 
result can be extended to 3D), Equation (7) can be expressed by two scalar equations: 
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Focusing on Equation (8), if the vertical acceleration a x and horizontal velocity v y as a function of y, can 
be estimated, then Equation (8) can be solved for the path P(y) by defining s = dP/dy = dxldy: 
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Equation (10) can be integrated once to solve for s(y ): 


V 


= \^dy' +Cl 

J v ( V ) 


( 11 ) 


Integrating once again yields the trajectory path function with two free boundary value constants, C\ and 

C 2 . 
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The constants of integration can be found by specifying the initial starting point and the slope of the initial 
starting point. Specify the slope and solve for C\i 


s(y 0 ) = s o 

Then solve for C 2 by specifying the initial starting value: 

P(y 0 ) = x o 


(14) 
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Application to CFD (Fluent™) Based Trajectories 

If the vertical acceleration a x as a function of horizontal distance y is approximated as a power law of the 
form: 

a x {y)~ccy- p - g , (16) 

where g is gravity, and while approximating the horizontal velocity v y as a constant v 0 , the double 
integral of Equation (13) can be evaluated. If the power law exponent (j is set equal to 2, then a further 
simplified form of Equation (12) results: 


Applying the boundary condition of Equation (14) results in 
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Finally, integrating Equation (13) gives the path P(y): 
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Applying the boundary condition of Equation (15) results in the following: 
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Combining Equations (18), (19), and (20) gives a complete trajectory path equation: 
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Comparison of Equation (21) with the Equation (A-l) describing a previous empirical model 
shows that most parameters are with the exception of a in this model and b in the empirical 
model. In the empirical model, the quantity bxjyo describes the tangent of the final trajectory 
angle in the limiting case of g = 0, which is approximately the case for large horizontal velocities 
within the CFD simulation domain: 


6 = tan 1 
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The parameter a in Equation (21) can be expressed in terms of b by taking the limit of the 
derivative of the path function as y goes to infinity with g = 0: 
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Solving for a and substituting into Equation (21) results in a trajectory path formula similar to the 
previous empirical formula, but now based on a true physical model: 


(24) 
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Results and Discussion 

The model described by Equation (24) is based on the assumption that the vertical acceleration o x as a 
function of horizontal distance y is approximated by a power law as described by Equation (16), along 
with a horizontal velocity that is described by a constant v 0 . For a practical application of this model, 
it is not necessary to find a fit to the vertical acceleration since the method previously used in the 
empirical model fitting can again be applied to this case. In this way, the particle trajectory path is 
determined by two primary fitting parameters, s 0 and b. For very small horizontal velocities, v 0 may also 
become a fitting parameter. 





Figure 3. Application of the parameter matrix interpolation of Equation (4) using the Imperical Trajectory 
Model of Equation (1), operating on the output of a Fluent CFD simulation of rocket with a nozzle height 
of 5 ft and a thrust of 67 N. Top: particle D = 1 jiim; bottom: particle D = 1000 jum. 




Figure 4. Solid lines: Full particle trajectories generated from post-processing of CFD simulations; dashed 
lines: trajectories generated by fitting the model of Equation (24) to solid lines. Top: D = 1 pm and x 0 = 

0.01 m; bottom: D = 1000 pm and x 0 = 0.01 m. 


With the example parameter interpolation matrix of Table 1, Equations (4) and (5) define a procedure to 
estimate the trajectory path of particles using either the empirical trajectory model of Equation (1) or the 
physical trajectory model of Equation (24). The number and extent of model variables is determined by 
the extent of the parameter interpolation matrix. In this example, only three particle variables were 
considered for simplicity of description, D, x 0 , and y 0 - Figure 3 shows the application of the parameter 


matrix interpolation of Equation (4) using the ETP model of Equation (1), operating on the output of a 
Fluent CFD simulation of a rocket with a nozzle height of 5 ft and a thrust of 67 N. The top set of plots 
corresponds to a particle of D = 1 pm, while the bottom plot is for a D - 1000 pm particle. The left plot 
is the parameter v y , corresponding to a constant horizontal velocity. The right plot is the gravity free 
trajectory angle, defined by Equation (2). The stair-stepped appearance of the plot is a consequnece of 
the coarse matrix used (Table 1) with N = 21. A larger matrix would smooth out the plot. Also, the 
simple Shepard’s interpolation defined by Equation (4) is known to generate non-smooth surfaces since 
the gradient at all reference points (all point defined by P,) is equal to zero because of Equation (4). A 
quadratic variation of Shepard’s interpolation (Renka, 1988) corrects this problem, but at the expense of a 
more complex algorithm. 

Conclusions 

In the process of verifying and characterizing the physical trajectory path model with Fluent test cases and 
comparing it to the previous empirical trajectory path model, the PTP model can be characterized by the 
following points: 

• The PTP model has a physical basis, starting with simplified estimates of vertical acceleration 
and horizontal velocity. The input parameters are x 0 and yo (initial starting point); s 0 derivative of 
trajectory dx(y)/dy, evaluated at y = y 0 ; and b where bxo/y 0 is the final trajectory angle before 
gravity pulls the particle down. Note that s 0 and b are the primary curve-fitting parameters. The 
final parameter is v 0 , an assumed constant horizontal velocity and is a secondary curve fitting 
parameter. 

• This model is time independent, showing vertical position rasa function of horizontal distance y, 
by x(y) = fly), as described by Equation (24). The first term on the right is due to simple ballistics 
and a spherically expanding rarefied gas so that the trajectory is a straight line intersecting (0,0,0), 
which is the point directly under the nozzle center on the surface. This term is generally the < 3° 
trajectory angle. The second term on the right is due to vertical acceleration, which may be 
positive or negative. The last term on the right is the gravity term, which for a particle with 
velocities less than escape velocity will eventually bring the particle back to the ground. 

• This model, even though it has the same number of parameters as the ETP model, seems to better 
describe the particle trajectories, even though it is far from perfect. Perfection in this context is 
traded for simplicity. 
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